﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "romb")]
    public static unsafe double romb(double a, double b, double eps, IntPtr f_x_ptr)
    {
        f_x = Marshal.GetDelegateForFunctionPointer<delegatefunc_x>(f_x_ptr);

        return romb(a, b, eps);
    }

    /// <summary>
    /// Romberg求积法
    /// f计算被积函数值f(x)的函数名。
    /// </summary>
    /// <param name="a">积分下限。</param>
    /// <param name="b">积分上限。要求b>a。</param>
    /// <param name="eps">积分精度要求。</param>
    /// <returns>函数返回积分值。</returns>
    public static unsafe double romb(double a, double b, double eps)
    {
        int m, n, i, k;
        double h, ep, p, x, s, q = 0.0;
        double[] y = new double[10];

        h = b - a;
        y[0] = h * (f_x(a) + f_x(b)) / 2.0;
        m = 1;
        n = 1;
        ep = eps + 1.0;
        while ((ep >= eps) && (m <= 9))
        {
            p = 0.0;
            for (i = 0; i <= n - 1; i++)
            {
                x = a + (i + 0.5) * h;
                p = p + f_x(x);
            }
            p = (y[0] + h * p) / 2.0;
            s = 1.0;
            for (k = 1; k <= m; k++)
            {
                s = 4.0 * s;
                q = (s * p - y[k - 1]) / (s - 1.0);
                y[k - 1] = p;
                p = q;
            }
            ep = Math.Abs(q - y[m - 1]);
            m = m + 1;
            y[m - 1] = q;
            n = n + n;
            h = h / 2.0;
        }
        return (q);
    }

    /*
    // Romberg求积法例
      int main()
      { 
          double a,b,eps,t,rombf(double);
          a=0.0; b=1.0; eps=0.000001;
          t = romb(a,b,eps,rombf);
          cout <<"t = " <<t <<endl;
          return 0;
      }
    // 计算被积函数值
      double rombf(double x)
      { 
          return(x/(4.0+x*x));
      }
    */
}

